# load necessary libraries
Warning message:
R graphics engine version 15 is not supported by this version of RStudio. The Plots tab will be disabled until a newer version of RStudio is installed.
library(Seurat)
The legacy packages maptools, rgdal, and rgeos, underpinning this package
will retire shortly. Please refer to R-spatial evolution reports on
https://r-spatial.org/r/2023/05/15/evolution4.html for details.
This package is now running under evolution status 0
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching SeuratObject
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(ggplot2)
library(CelltypeR)
Warning: replacing previous import ‘data.table::last’ by ‘dplyr::last’ when loading ‘CelltypeR’
Warning: replacing previous import ‘data.table::first’ by ‘dplyr::first’ when loading ‘CelltypeR’
Warning: replacing previous import ‘data.table::between’ by ‘dplyr::between’ when loading ‘CelltypeR’
Warning: replacing previous import ‘dplyr::filter’ by ‘flowCore::filter’ when loading ‘CelltypeR’
Warning: replacing previous import ‘ggplot2::margin’ by ‘randomForest::margin’ when loading ‘CelltypeR’
Warning: replacing previous import ‘dplyr::combine’ by ‘randomForest::combine’ when loading ‘CelltypeR’
Warning: replacing previous import ‘flowCore::filter’ by ‘dplyr::filter’ when loading ‘CelltypeR’
Warning: replacing previous import ‘flowViz::contour’ by ‘graphics::contour’ when loading ‘flowStats’
Attaching package: ‘CelltypeR’
The following object is masked from ‘package:ggplot2’:
annotate
library(flowCore)
Read in the flow data This data should be the gated live cells.
All samples need to be in one folder. Here we have
sampleNames(flowset)
[1] "MBO_154d_AIW_n1_Live_Live.fcs"
[2] "MBO_154d_AIW_n2_Live_Live.fcs"
[3] "MBO_154d_AIW_n3_Live_Live.fcs"
[4] "MBO_154d_AIW_n4_Live_Live.fcs"
[5] "MBO_30d_AIW_n1_Single Cells_Live_Live.fcs"
[6] "MBO_30d_AIW_n2_Single Cells_Live_Live.fcs"
[7] "MBO_30d_AIW_n3_Single Cells_Live_Live.fcs"
[8] "MBO_30d_AIW_n4_Single Cells_Live_Live.fcs"
[9] "MBO_60d_AIW_n1_Live_Live.fcs"
[10] "MBO_60d_AIW_n2_Live_Live.fcs"
[11] "MBO_60d_AIW_n3_Live_Live.fcs"
[12] "MBO_60d_AIW_n4_Live_Live.fcs"
[13] "MBO_97d_AIW_n1_Live_Live.fcs"
[14] "MBO_97d_AIW_n2_Live_Live.fcs"
[15] "MBO_97d_AIW_n3_Live_Live.fcs"
[16] "MBO_97d_AIW_n4_Live_Live.fcs"
Rename samples with shorter names
sampleNames(flowset) <- sampleNames(flowset) <- c("AIW_150_R1","AIW_150_R2","AIW_150_R3","AIW_150_R4",
"AIW_30_R1","AIW_30_R2","AIW_30_R3","AIW_30_R4",
"AIW_60_R1","AIW_60_R2","AIW_60_R3","AIW_60_R4",
"AIW_100_R1","AIW_100_R2","AIW_100_R3","AIW_100_R4"
)
sampleNames(flowset)
[1] "AIW_150_R1" "AIW_150_R2" "AIW_150_R3" "AIW_150_R4"
[5] "AIW_30_R1" "AIW_30_R2" "AIW_30_R3" "AIW_30_R4"
[9] "AIW_60_R1" "AIW_60_R2" "AIW_60_R3" "AIW_60_R4"
[13] "AIW_100_R1" "AIW_100_R2" "AIW_100_R3" "AIW_100_R4"
Harmonize data to remove batch or technical variation
This requires us to look and see where there are two peaks to align. We need to visualize the peaks of the transformed data before aligning.
colnames(df)
[1] "FSC" "FSC" "FSC" "SSC" "SSC"
[6] "SSC" "CD56" "CD56" "CD56" "SSEA4"
[11] "SSEA4" "SSEA4" "CD140a" "CD140a" "CD140a"
[16] "CD29" "CD29" "CD29" "CD44" "CD44"
[21] "CD44" "LiveDead" "LiveDead" "LiveDead" "CD184"
[26] "CD184" "CD184" "CD24" "CD24" "CD24"
[31] "CD15" "CD15" "CD15" "TH" "TH"
[36] "TH" "CD49f" "CD49f" "CD49f" "CD133"
[41] "CD133" "CD133" "Time" "Sample" "cell"
Now we have made all the different processing of the fsc files. We can visualize the intensity in cell density plots to see the alignment
#plotdensity_flowset(flowset)
plotdensity_flowset(flowset_biexp)
Warning in melt(lapply(as.list(flowset@frames), function(x) { :
The melt generic in data.table has been passed a list and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(lapply(as.list(flowset@frames), function(x) { x = as.data.frame(x@exprs)})). In the next version, this warning will become an error.
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
Warning in geom_density_ridges(alpha = 0.4, verbose = FALSE) :
Ignoring unknown parameters: `verbose`
Picking joint bandwidth of 0.106
Picking joint bandwidth of 0.0931
Picking joint bandwidth of 0.0233
Picking joint bandwidth of 0.124
Picking joint bandwidth of 0.107
Picking joint bandwidth of 0.023
Picking joint bandwidth of 0.224
Picking joint bandwidth of 0.202
Picking joint bandwidth of 0.396
Picking joint bandwidth of 0.236
Picking joint bandwidth of 0.166
Picking joint bandwidth of 0.199
Picking joint bandwidth of 0.25
Picking joint bandwidth of 0.232
Picking joint bandwidth of 0.374
Picking joint bandwidth of 0.227
Picking joint bandwidth of 0.212
Picking joint bandwidth of 0.352
Picking joint bandwidth of 0.281
Picking joint bandwidth of 0.24
Picking joint bandwidth of 0.325
Picking joint bandwidth of 0.213
Picking joint bandwidth of 0.152
Picking joint bandwidth of 0.172
Picking joint bandwidth of 0.351
Picking joint bandwidth of 0.265
Picking joint bandwidth of 0.407
Picking joint bandwidth of 0.246
Picking joint bandwidth of 0.236
Picking joint bandwidth of 0.195
Picking joint bandwidth of 0.234
Picking joint bandwidth of 0.222
Picking joint bandwidth of 0.394
Picking joint bandwidth of 0.309
Picking joint bandwidth of 0.281
Picking joint bandwidth of 0.398
Picking joint bandwidth of 0.343
Picking joint bandwidth of 0.275
Picking joint bandwidth of 0.337
Picking joint bandwidth of 0.236
Picking joint bandwidth of 0.205
Picking joint bandwidth of 0.34
Picking joint bandwidth of 0.14
plotdensity_flowset(flowset_align)
Warning in melt(lapply(as.list(flowset@frames), function(x) { :
The melt generic in data.table has been passed a list and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(lapply(as.list(flowset@frames), function(x) { x = as.data.frame(x@exprs)})). In the next version, this warning will become an error.
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
No id variables; using all as measure variables
Warning in geom_density_ridges(alpha = 0.4, verbose = FALSE) :
Ignoring unknown parameters: `verbose`
Picking joint bandwidth of 0.105
Picking joint bandwidth of 0.0921
Picking joint bandwidth of 0.0231
Picking joint bandwidth of 0.123
Picking joint bandwidth of 0.106
Picking joint bandwidth of 0.023
Picking joint bandwidth of 0.224
Picking joint bandwidth of 0.202
Picking joint bandwidth of 0.395
Picking joint bandwidth of 0.232
Picking joint bandwidth of 0.164
Picking joint bandwidth of 0.201
Picking joint bandwidth of 0.251
Picking joint bandwidth of 0.232
Picking joint bandwidth of 0.372
Picking joint bandwidth of 0.226
Picking joint bandwidth of 0.208
Picking joint bandwidth of 0.35
Picking joint bandwidth of 0.269
Picking joint bandwidth of 0.238
Picking joint bandwidth of 0.325
Picking joint bandwidth of 0.213
Picking joint bandwidth of 0.152
Picking joint bandwidth of 0.172
Picking joint bandwidth of 0.351
Picking joint bandwidth of 0.266
Picking joint bandwidth of 0.409
Picking joint bandwidth of 0.247
Picking joint bandwidth of 0.238
Picking joint bandwidth of 0.196
Picking joint bandwidth of 0.235
Picking joint bandwidth of 0.223
Picking joint bandwidth of 0.393
Picking joint bandwidth of 0.303
Picking joint bandwidth of 0.281
Picking joint bandwidth of 0.399
Picking joint bandwidth of 0.338
Picking joint bandwidth of 0.272
Picking joint bandwidth of 0.336
Picking joint bandwidth of 0.232
Picking joint bandwidth of 0.202
Picking joint bandwidth of 0.337
Picking joint bandwidth of 0.14
#plotdensity_flowset(flowset_retro)
AB <- c("TH","CD24","CD56","CD29","CD15","CD184","CD133","SSEA4","CD44","CD49f","CD140a")
df_select <- df %>% select(c("TH","CD24","CD56","CD29","CD15","CD184","CD133","SSEA4","CD44","CD49f","CD140a","Sample"))
seu <- make_seu(df_select, AB_vector = AB)
Centering and scaling data matrix
|
| | 0%
|
|====================================================================| 100%
Warning in irlba(A = t(x = object), nv = npcs, ...) :
You're computing too large a percentage of total singular values, use a standard svd instead.
Warning: Requested number is larger than the number of available items (11). Setting to 11.
Warning: Requested number is larger than the number of available items (11). Setting to 11.
Warning: Requested number is larger than the number of available items (11). Setting to 11.
Warning: Requested number is larger than the number of available items (11). Setting to 11.
Warning: Requested number is larger than the number of available items (11). Setting to 11.
PC_ 1
Positive: CD24, CD15, CD29, CD140a, CD56, CD133
Negative: TH, SSEA4, CD49f, CD44, CD184
PC_ 2
Positive: CD15, CD24, CD140a, CD56, TH, SSEA4
Negative: CD44, CD133, CD29, CD49f, CD184
PC_ 3
Positive: CD44, CD15, CD24, CD140a, CD29, CD56
Negative: CD184, CD49f, CD133, TH, SSEA4
PC_ 4
Positive: CD184, CD44, CD15, CD24, TH, CD56
Negative: CD29, CD133, SSEA4, CD140a, CD49f
PC_ 5
Positive: CD56, CD49f, TH, CD133, CD24, CD15
Negative: SSEA4, CD184, CD140a, CD29, CD44
You can save the data frame and seurat object for later
# to save the df for later
write.csv(df, paste(output_path,"df2000Timecourse4.csv", sep = ""))
# save the seurat object
saveRDS(seu, paste(output_path, "seuratObjectTimecourseAIW_4.RDS", sep = ""))
Check out the data object
RidgePlot(seu, features = AB, log = TRUE)
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0762
Warning: Removed 709 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.061
Warning: Removed 311 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0553
Warning: Removed 235 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0564
Warning: Removed 230 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0572
Warning: Removed 488 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0872
Warning: Removed 1453 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0602
Warning: Removed 311 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0597
Warning: Removed 2300 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0621
Warning: Removed 713 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.0817
Warning: Removed 1088 rows containing non-finite values (`stat_density_ridges()`).
Warning in self$trans$transform(x) : NaNs produced
Warning: Transformation introduced infinite values in continuous x-axis
Picking joint bandwidth of 0.06
Warning: Removed 718 rows containing non-finite values (`stat_density_ridges()`).
DimPlot(seu, group.by = "Sample")
# note the dimensions (number of PC to include needs to be one less than the number in the antibody panel)
explore_param(input = seu,
cluster_method = "louvain",
df_input = df.input,
flow_k = NULL,
pheno_lou_kn = c(80,120,200),
lou_resolution = c(0.6),
pcdim = 1:10,
run.plot = TRUE,
run.stats = FALSE,
save_to = NULL)
Computing nearest neighbor graph
Computing SNN
13:40:53 UMAP embedding parameters a = 0.9922 b = 1.112
13:40:53 Read 32000 rows and found 10 numeric columns
13:40:53 Using Annoy for neighbor search, n_neighbors = 80
13:40:53 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:40:57 Writing NN index file to temp file /var/folders/k4/khtkczkd5tn732ftjpwgtr240000gn/T//RtmpcsjJYG/file2f9173d2fe46
13:40:57 Searching Annoy index using 1 thread, search_k = 8000
13:41:33 Annoy recall = 100%
13:41:34 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 80
13:41:39 Initializing from normalized Laplacian + noise (using irlba)
13:41:41 Commencing optimization for 200 epochs, with 3116580 positive edges
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:42:24 Optimization finished
[1] "finding neighbours for kn 80"
Warning: The following arguments are not used: reduction
Warning: The following arguments are not used: reduction
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 32000
Number of edges: 4000185
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8744
Number of communities: 13
Elapsed time: 26 seconds
[1] "completed louvain kn 80 resolution 0.6"
[1] 32000
Computing nearest neighbor graph
Computing SNN
13:43:41 UMAP embedding parameters a = 0.9922 b = 1.112
13:43:41 Read 32000 rows and found 10 numeric columns
13:43:41 Using Annoy for neighbor search, n_neighbors = 120
13:43:41 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:43:45 Writing NN index file to temp file /var/folders/k4/khtkczkd5tn732ftjpwgtr240000gn/T//RtmpcsjJYG/file2f912133e112
13:43:45 Searching Annoy index using 1 thread, search_k = 12000
13:44:34 Annoy recall = 100%
13:44:35 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 120
13:44:42 Initializing from normalized Laplacian + noise (using irlba)
13:44:46 Commencing optimization for 200 epochs, with 3950280 positive edges
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:45:33 Optimization finished
[1] "finding neighbours for kn 120"
Warning: The following arguments are not used: reduction
Warning: The following arguments are not used: reduction
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 32000
Number of edges: 5858818
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8655
Number of communities: 13
Elapsed time: 41 seconds
[1] "completed louvain kn 120 resolution 0.6"
[1] 32000
Computing nearest neighbor graph
Computing SNN
13:47:25 UMAP embedding parameters a = 0.9922 b = 1.112
13:47:25 Read 32000 rows and found 10 numeric columns
13:47:25 Using Annoy for neighbor search, n_neighbors = 200
13:47:25 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:47:29 Writing NN index file to temp file /var/folders/k4/khtkczkd5tn732ftjpwgtr240000gn/T//RtmpcsjJYG/file2f9154fd75c5
13:47:30 Searching Annoy index using 1 thread, search_k = 20000
13:48:48 Annoy recall = 100%
13:48:49 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 200
13:49:01 Initializing from normalized Laplacian + noise (using irlba)
13:49:06 Commencing optimization for 200 epochs, with 4802826 positive edges
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:49:59 Optimization finished
[1] "finding neighbours for kn 200"
Warning: The following arguments are not used: reduction
Warning: The following arguments are not used: reduction
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 32000
Number of edges: 9488963
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8514
Number of communities: 10
Elapsed time: 68 seconds
[1] "completed louvain kn 200 resolution 0.6"
[1] 32000
$cluster
NULL
# this function does not add the clustering to the seurat object it is to explore and run stats. If you choose save = output_path then you will save a seurat object with the clustering, one for each kn.
After checking parameters run cluster function to make add the clustering information into the
seu <- get_clusters(seu, method = "louvain",
df_input = df.input,
k = 80,
resolution = c(1),
plots = TRUE,
save_plots = FALSE)
Computing nearest neighbor graph
Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 32000
Number of edges: 4000185
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8450
Number of communities: 18
Elapsed time: 26 seconds
[1] "method is Louvain"
15:38:49 UMAP embedding parameters a = 0.4502 b = 1.076
15:38:49 Read 32000 rows and found 10 numeric columns
15:38:49 Using Annoy for neighbor search, n_neighbors = 80
15:38:49 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:38:54 Writing NN index file to temp file /var/folders/k4/khtkczkd5tn732ftjpwgtr240000gn/T//RtmpcsjJYG/file2f9146a08787
15:38:54 Searching Annoy index using 1 thread, search_k = 8000
15:39:30 Annoy recall = 100%
15:39:31 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 80
15:39:36 Initializing from normalized Laplacian + noise (using irlba)
15:39:39 Commencing optimization for 200 epochs, with 3116580 positive edges
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:40:21 Optimization finished
# if save_plots is TRUE the feature plots will be save as 13 plots in one file. The UMAP with cluster numbers will also be saved.
Visualize expression on UMAP and with heat maps
AB <- c("TH","CD24","CD56","CD29","CD15","CD184","CD133","SSEA4","CD44","CD49f","CD140a")
# the cluster labels will match the active ident
Idents(seu) <- "RNA_snn_res.1"
# this will let us see one at at time
for (i in AB) {
print(FeaturePlot(seu, features = i, min.cutoff = 'q1', max.cutoff = 'q97', label = TRUE))
}